install.packages("car")
install.packages("gmodels")
install.packages("xtable")
install.packages("effects")
install.packages("separationplot")
install.packages("pscl")
install.packages("MCMCpack")
install.packages("UsingR")
 #v.3.5
install.packages("gdata")
install.packages("gplots")
install.packages("foreign")
install.packages("lme4")
install.packages("nlme")
install.packages("texreg")
install.packages("stargazer")
install.packages("ggplot2")
install.packages("reshape2")
install.packages("ggcorrplot")
install.packages("data.table")



library(verification)
library(car)
library(gmodels)
library(xtable)
library(effects)
library(separationplot)
library(pscl)
library(MCMCpack)
library(UsingR)
library(Zelig) #v.3.5
library(gdata)
library(gplots)
library(foreign)
library(lme4)
library(nlme)
library(texreg)
library(stargazer)
library(ggplot2)
library(reshape2)
library(ggcorrplot)
library(data.table)

######
#read data
kos<-read.dta("C:\\Users\\Daniel\\Dropbox\\Journal (II)\\Replication Files\\Mirilovic\\Kosovo Dataset.dta")

######
#Table 1, column 1
l1<-glm(kosovo_recognition_dummy~NATO2007+socialism+gdppc2007+
          usfdigdp19901999_modified+dem2007+percent_muslim_arda+
          rr_additive_index+ Griffiths+ ethnic_diversity_alesina+
          belgrad_distance+ percent_albanian, family=binomial(link="logit"),
        data=kos)
summary(l1)


######
#Table 1, column 2
m1<-lmer(kosovo_recognition_dummy~NATO2007+socialism+gdppc2007+
           usfdigdp19901999_modified+dem2007+percent_muslim_arda+
           rr_additive_index+ Griffiths+ ethnic_diversity_alesina+
           belgrad_distance+ percent_albanian+ (1|regnum), REML=T, data=kos)
summary(m1)
summary(kos$rr_additive_index)
0.02203/(0.02203+0.17492)


##calculate ICC
m2<-lmer(kosovo_recognition_dummy~ (1|regnum), REML=T, data=kos)
summary(m2)
0.05294/(0.05294+0.21206)

kos$NATO01[kos$NATO2007=="No"]<-0
kos$NATO01[kos$NATO2007=="Yes"]<-1

######
#Table 1, column 3 (in stata)
xtmelogit kosovo_recognition_dummy rr_additive_index percent_muslim_arda NATO2007 socialism gdppc2007  usfdigdp19901999_modified dem2007 
Griffiths ethnic_diversity_alesina belgrad_distance percent_albanian || regnum:, var mle


######
###Marginal effects based on model l1
#religious regulation
dat <- data.frame(rr_additive_index = rep(seq(from = min(kos$rr_additive_index,na.rm=T), 
                                              to = max(kos$rr_additive_index,na.rm=T), length.out = 1000)),
                  NATO01 =mean(kos$NATO01),
                  socialism =mean(kos$socialism),
                  gdppc2007 =mean(kos$gdppc2007,na.rm=T),
                  usfdigdp19901999_modified =mean(kos$usfdigdp19901999_modified,na.rm=T),
                  dem2007=mean(kos$dem2007,na.rm=T),
                  Griffiths=mean(kos$Griffiths,na.rm=T),
                  ethnic_diversity_alesina=mean(kos$ethnic_diversity_alesina,na.rm=T),
                  belgrad_distance=mean(kos$belgrad_distance,na.rm=T),
                  percent_muslim_arda=mean(kos$percent_muslim_arda,na.rm=T),
                  percent_albanian=mean(kos$percent_albanian,na.rm=T))


ndat <- cbind(dat, predict(l1, newdata=dat, type = "response",se.fit=T)) 
newdat <- melt(ndat, id.vars = c("rr_additive_index","fit","se.fit"))
newdat2 <- newdat[1:1000,]
newdat2$ucl <- newdat2$fit + 1.96 * newdat2$se.fit
newdat2$lcl <- newdat2$fit - 1.96 * newdat2$se.fit


tail(newdat2)                                    
b<-ggplot(newdat2, aes(x = rr_additive_index, y = fit)) +
  geom_smooth(aes(ymin = lcl, ymax = ucl), stat="identity")+
  labs(title = "", x="Religious Regulation", 
       y="Probability of Recognition")
theme(legend.position="bottom")+
  scale_x_continuous(breaks=c(2.0,2.1,2.2,2.3,2.4,2.5), 
                     labels=c("0", "10", "20", "30", "40","50"))
b
ggsave(file="relregulation.pdf",width=7,height=7)

#religious ties/percent muslim
dat <- data.frame(percent_muslim_arda = rep(seq(from = min(kos$percent_muslim_arda,na.rm=T), 
                                                to = max(kos$percent_muslim_arda,na.rm=T), length.out = 1000)),
                  NATO01 =mean(kos$NATO01),
                  socialism =mean(kos$socialism),
                  gdppc2007 =mean(kos$gdppc2007,na.rm=T),
                  usfdigdp19901999_modified =mean(kos$usfdigdp19901999_modified,na.rm=T),
                  dem2007=mean(kos$dem2007,na.rm=T),
                  Griffiths=mean(kos$Griffiths,na.rm=T),
                  ethnic_diversity_alesina=mean(kos$ethnic_diversity_alesina,na.rm=T),
                  belgrad_distance=mean(kos$belgrad_distance,na.rm=T),
                  rr_additive_index=mean(kos$rr_additive_index,na.rm=T),
                  percent_albanian=mean(kos$percent_albanian,na.rm=T))

ndat <- cbind(dat, predict(l1, newdata=dat, type = "response",se.fit=T)) 
newdat <- melt(ndat, id.vars = c("percent_muslim_arda","fit","se.fit"))
newdat2 <- newdat[1:1000,]
newdat2$ucl <- newdat2$fit + 1.96 * newdat2$se.fit
newdat2$lcl <- newdat2$fit - 1.96 * newdat2$se.fit


tail(newdat2)                                    
b<-ggplot(newdat2, aes(x = percent_muslim_arda, y = fit)) +
  geom_smooth(aes(ymin = lcl, ymax = ucl), stat="identity")+
  labs(title = "", x="Percent Muslim", 
       y="Probability of Recognition")
theme(legend.position="bottom")+
  scale_x_continuous(breaks=c(2.0,2.1,2.2,2.3,2.4,2.5), 
                     labels=c("0", "10", "20", "30", "40","50"))
b
ggsave(file="percentmuslim.pdf",width=7,height=7)


######
#Appendix Tables and Figures

######
#Correlation matrix - appendix figure 1
myvars<-c("NATO01","socialism","gdppc2007",
          "usfdigdp19901999_modified","dem2007","percent_muslim_arda",
          "rr_additive_index", "Griffiths", "ethnic_diversity_alesina",
          "belgrad_distance", "percent_albanian")

kos2<-kos[myvars]
corr <- round(cor(kos2), 1)
corr
xtable(corr)
setnames(kos2,"NATO01","NATO")
setnames(kos2,"socialism","Socialism")
setnames(kos2,"gdppc2007","GDPpc")
setnames(kos2,"usfdigdp19901999_modified","USFDI")
setnames(kos2,"dem2007","Democracy")
setnames(kos2,"percent_muslim_arda","Muslim")
setnames(kos2,"rr_additive_index","Regulation")
setnames(kos2,"Griffiths","Vulnerable")
setnames(kos2,"ethnic_diversity_alesina","Diversity")
setnames(kos2,"belgrad_distance","Distance")
setnames(kos2,"percent_albanian","Albanian")

ggcorrplot(corr, hc.order = TRUE,
           lab = F,legend.title = "")
ggsave("corr_nolab.pdf",width=7,height=7)
dev.off()

######
#NATO membership effect - appendix fig 2a

lo<-setx(z1, NATO2007 = "No")
hi<-setx(z1, NATO2007 = "Yes")

slo<- sim(z1, x = lo)
shi<- sim(z1, x = hi)

summary(slo)
summary(shi)

low<-melt(slo$qi$ev)
high<-melt(shi$qi$ev)

ll<-data.frame(level=rep("No NATO Membership"), low)
hh<-data.frame(level=rep("NATO Membership"), high)

hl<-rbind(hh,ll)

g<-ggplot(hl, aes(x=level, y=value))
g+  geom_boxplot()+
  ylab("Probability of Recognition") +
  xlab("")+
  ylim(0,1)
ggsave("nato01.pdf",width=7,height=7)
dev.off()

######
##domestic vulnerability effect  - appendix fig 2b

lo<-setx(z1, Griffiths = 0)
hi<-setx(z1, Griffiths = 1)

slo<- sim(z1, x = lo)
shi<- sim(z1, x = hi)

low<-melt(slo$qi$ev)
high<-melt(shi$qi$ev)

ll<-data.frame(level=rep("Not Vulnerable"), low)
hh<-data.frame(level=rep("Vulnerable"), high)

hl<-rbind(hh,ll)

g<-ggplot(hl, aes(x=level, y=value))
g+  geom_boxplot()+
  ylab("Probability of Recognition") +
  xlab("")+
  ylim(0,1)
ggsave("vulnerable01.pdf",width=7,height=7)
dev.off()

######
#vif test - appendix table 1

vif(l1) 

######
#Anova test - appendix table 2
l2<-glm(kosovo_recognition_dummy~NATO2007+socialism+gdppc2007+
          usfdigdp19901999_modified+dem2007+
          Griffiths+ ethnic_diversity_alesina+
          belgrad_distance+ percent_albanian, family=binomial(link="logit"),
        data=kos)
summary(l2)
a<-anova(m1,test = "Chisq")
xtable(a)
texreg(a)
texreg(list(l1,m1),booktabs = TRUE, dcolumn = TRUE)

z1<-zelig(kosovo_recognition_dummy~NATO2007+socialism+gdppc2007+
            usfdigdp19901999_modified+dem2007+percent_muslim_arda+
            rr_additive_index+ Griffiths+ ethnic_diversity_alesina+
            belgrad_distance+ percent_albanian, model="logit",
          data=kos)
summary(z1)
